def determine_defect(fp,log_fp):
    import GlobalParam
    import sys
    gp=GlobalParam.iparam

    # Neighbor cells
    # Time step
    eatline=fp.readline()
    eatline=fp.readline()

    # Number of atoms
    eatline=fp.readline()
    eatline=fp.readline()
    natom=int(eatline)

    if natom != gp.natom:
        print("File (%s) are not consistent with initial source file." % src_file_name)
        raise Exception

    # Read dummy lines
    for i in range(5):
        eatline=fp.readline()

    large=1.0e+5    

    # Initialize all counters
    gp.reset_count()

    showprocess=0.0
    processdelta=0.1
    # Iteration for each atom
    log_str=" -complete reading %3.2f %%\n"
    for i in range(gp.get_natom()):
        curprocess=float(i)/float(gp.get_natom())
        if curprocess>=showprocess:
            log_fp.write(log_str % (showprocess*100))
            sys.stdout.write(log_str % (showprocess*100))
            showprocess+=processdelta
        eatline=fp.readline()
        # output format: atom_id atom_type px py pz
        eatword=eatline.split()
        atom_id=int(eatword[0])
        atom_type=int(eatword[1])
        px=float(eatword[2])
        py=float(eatword[3])
        pz=float(eatword[4])
        # Check if x position is out of box
        if px<0:
            px=px+gp.get_boxwidth()
        elif px>=gp.get_boxwidth():
            px=px-gp.get_boxwidth()
        else:
            pass
        # Check if y position is out of box
        if py<0:
            py=py+gp.get_boxwidth()
        elif py>=gp.get_boxwidth():
            py=py-gp.get_boxwidth()
        else:
            pass
        # Check if z position is out of box
        if pz<0:
            pz=pz+gp.get_boxwidth()
        elif pz>=gp.get_boxwidth():
            pz=pz-gp.get_boxwidth()
        else:
            pass

        dist=large

        # Find which cell does it belong?
        cellx=int(px/gp.get_lattice_param())
        celly=int(py/gp.get_lattice_param())
        cellz=int(pz/gp.get_lattice_param())

        # Find nearest lattice point
        basis_no=-1
        for item1 in gp.neighboriter(cellx,celly,cellz):
            for item2 in gp.gen_lattice_point_in_cell(item1):
                cdist=(item2[0]-px)**2+\
                      (item2[1]-py)**2+\
                      (item2[2]-pz)**2
                if dist>cdist:
                    # save at nearest lattice
                    ix=item2[3]
                    iy=item2[4]
                    iz=item2[5]
                    ppx=px
                    ppy=py
                    ppz=pz
                    basis_no=item2[6]
                    dist=cdist
        if basis_no<0:
            print("Fatal error: Can't find nearest lattice point.")
            raise Exception
        else:
            gp.count_nearest_point(ix,iy,iz,ppx,ppy,ppz,basis_no)

def index2xyz(n):
    import GlobalParam
    gp=GlobalParam.iparam
    iz=n/gp.get_ncell_in_plane()
    iy=(n-iz*gp.get_ncell_in_plane())/gp.get_ncell_in_line()
    ix=n-iz*gp.get_ncell_in_plane()-iy*gp.get_ncell_in_line()
    x0=float(ix)*gp.get_lattice_param()
    y0=float(iy)*gp.get_lattice_param()
    z0=float(iz)*gp.get_lattice_param()
    return ix,iy,iz,x0,y0,z0

def count_defect():
    import GlobalParam
    gp=GlobalParam.iparam
    for i in range(gp.get_ncell()):
        ix,iy,iz,x0,y0,z0=index2xyz(i)
        for j in range(gp.get_nbasis()):
            x=x0+gp.get_basis2_pos(j)[0]
            y=y0+gp.get_basis2_pos(j)[1]
            z=z0+gp.get_basis2_pos(j)[2]
            # natom=gp.unitcell[i].get_basepoint(j)
            natom=gp.unitcell[i].get_natom(j)
            if natom == 0:
                gp.count_vacancy()
                gp.add_vacancy_list(i,j)
                gp.n_defect_atom+=1
            elif natom == 1:
                gp.count_normal()
            elif natom == 2:
                gp.count_SIA()
                gp.add_SIA_list(i,j)
                gp.n_defect_atom+=2
            else:
                gp.count_odd()
                gp.add_odd_list(i,j)
                gp.n_defect_atom+=natom

def write2log(fp,i):
    import GlobalParam
    gp=GlobalParam.iparam
    fp.write("%8d %5d %5d %5d %5d\n" % \
      (gp.file_number_list[i],gp.get_nvacancy(),gp.get_nSIA(),gp.get_nodd(),gp.get_nnormal()))
    print("  File[%3d].. name: %15s V=%5d  SIA=%5d  Odd=%5d  Normal=%5d" % \
      (i,gp.get_srcfile(i),gp.get_nvacancy(),gp.get_nSIA(),gp.get_nodd(),gp.get_nnormal()))


def write2out(fp,i):
    import GlobalParam
    gp=GlobalParam.iparam
    fp.write("%8d %5d %5d %5d %5d\n" % \
      (gp.file_number_list[i],gp.get_nvacancy(),gp.get_nSIA(),gp.get_nodd(),gp.get_nnormal()))

def write2xml(fp,str1="C",str2="I",str3="Fe"):
    import GlobalParam
    gp=GlobalParam.iparam
    # Number of atoms
    fp.write("%d\n" % (gp.get_n_defect_atom()+gp.get_nxmlframe()))
    # Comment
    fp.write("US Fe\n")
    # Write frames
    for item in gp.gen_xmlframe():
        fp.write("%s %f %f %f\n" % (str1,item[0],item[1],item[2]))
    # Write defects
    for i in range(gp.get_nvacancy()):
        cell_index=gp.vacancy_list[i][0]
        base_index=gp.vacancy_list[i][1]
        ix,iy,iz,x0,y0,z0=index2xyz(cell_index)
        x=x0+gp.get_basis2_pos(base_index)[0]
        y=y0+gp.get_basis2_pos(base_index)[1]        
        z=z0+gp.get_basis2_pos(base_index)[2]
        fp.write("%s %f %f %f\n" % (str2,x,y,z))
    for i in range(gp.get_nSIA()):
        cell_index=gp.vacancy_list[i][0]
        base_index=gp.vacancy_list[i][1]
        for j in range(gp.unitcell[cell_index].get_natom(base_index)):
            x,y,z=gp.unitcell[cell_index].get_pos(base_index,j)
            fp.write("%s %f %f %f"% (str3,x,y,z))
    for i in range(gp.get_nodd()):
        cell_index=gp.odd_list[i][0]
        base_index=gp.odd_list[i][1]
        for j in range(gp.unitcell[cell_index].get_natom(base_index)):
            x,y,z=gp.unitcell[cell_index].get_pos(base_index,j)
            fp.write("%s %f %f %f"% (str3,x,y,z))

def write2kmc(fp):
    import GlobalParam
    gp=GlobalParam.iparam
    str1="PKA_E %s eV %d\n"
    # Write defects
    fp.write(str1 % (gp.get_case_string(),gp.get_nvacancy()+gp.get_nSIA()+gp.get_nodd()))
    dcount=1
    for i in range(gp.get_nvacancy()):
        cell_index=gp.vacancy_list[i][0]
        base_index=gp.vacancy_list[i][1]
        ix,iy,iz,x0,y0,z0=index2xyz(cell_index)
        x=x0+gp.get_basis2_pos(base_index)[0]
        y=y0+gp.get_basis2_pos(base_index)[1]        
        z=z0+gp.get_basis2_pos(base_index)[2]
        fp.write("%d V1 %f %f %f\n" % (dcount,x,y,z))
        dcount+=1
    for i in range(gp.get_nSIA()):
        cell_index=gp.SIA_list[i][0]
        base_index=gp.SIA_list[i][1]
        xx=[]
        yy=[]
        zz=[]
        for j in range(gp.unitcell[cell_index].get_natom(base_index)):
            x,y,z=gp.unitcell[cell_index].get_pos(base_index,j)
            xx.append(x)
            yy.append(y)
            zz.append(z)
        dx=xx[0]-xx[1]
        dy=yy[0]-yy[1]
        dz=zz[0]-zz[1]
        if dx>0:
            d1=1
        elif dx<0:
            d1=-1
        else:
            d1=0
        if dy>0:
            d2=1
        elif dy<0:
            d2=-1
        else:
            d2=0
        if dz>0:
            d3=1
        elif dz<0:
            d3=-1
        else:
            dz=0
        fp.write("%d I1 %f %f %f %d %d %d\n" % (dcount,x,y,z,d1,d2,d3))
        dcount+=1
    for i in range(gp.get_nodd()):
        cell_index=gp.odd_list[i][0]
        base_index=gp.odd_list[i][1]
        ix,iy,iz,x0,y0,z0=index2xyz(cell_index)
        x=x0+gp.get_basis2_pos(base_index)[0]
        y=y0+gp.get_basis2_pos(base_index)[1]        
        z=z0+gp.get_basis2_pos(base_index)[2]
        fp.write("%d SS %f %f %f\n" % (dcount,x,y,z))
        dcount+=1

def countbysource(out_fp,log_fp):
    import GlobalParam
    import sys
    gp=GlobalParam.iparam
    log_fp.write("  - Number of defects for each file..\n")
    log_fp.write("Time     Vacancy  SIA  Odd Normal\n")
    out_fp.write("  - Number of defects for each file..\n")
    out_fp.write("Time     Vacancy  SIA  Odd Normal\n")
    str1=" Reading source file... [%d/%d] %s\n"
    for i in range(gp.get_nsrcfile()):
        # Default
        source_fp=open(gp.get_srcfile(i),'r')
        sys.stdout.write(str1 % (i,gp.get_nsrcfile(),gp.get_srcfile(i)))
        log_fp.write(str1 % (i,gp.get_nsrcfile(),gp.get_srcfile(i)))
        determine_defect(source_fp,log_fp)
        count_defect()
        write2log(log_fp,i)
        write2out(out_fp,i)
        source_fp.close()

        # Additional writing
        xml_fp=open(gp.get_xmlfile(i),'w')
        write2xml(xml_fp)
        xml_fp.close()

        kmc_fp=open(gp.get_kmcfile(i),'w')
        write2kmc(kmc_fp)
        kmc_fp.close()
